;  calculating matrix making
;  rotary around arbitrary vector
;  implemented by Maciej Guba,
;  macgub3@wp.pl , http://macgub.co.pl
;  ebx - normalized vector 3 dwords float
;  esi - cosinus and sinus of angle
;  edi - desired matrix


make_arbitrary_mx:
 .nx	  equ [ebx]
 .ny	  equ [ebx+4]
 .nz	  equ dword[ebx+8]
 .sinus   equ [esi]
 .cosinus equ dword[esi+4]

 or	  eax,-1
 push	  eax
 neg	  eax
 push	  eax
 cvtpi2ps xmm0,[esp]   ; xm0 = lo->hi 1.0 -1.0
 movups   xmm1,.nx
 movaps   xmm5,xmm1    ; xm5 = nxyz?
 movaps   xmm6,xmm1
 mulps	  xmm1,xmm1    ; xm1 = nxyz?^2
 movaps   xmm2,xmm0
 shufps   xmm2,xmm2,0  ; xm2 = brdcst 1.0
 movaps   xmm4,xmm2
 subps	  xmm2,xmm1    ; xm2 = 1 - nxyz?^2

 movlps   xmm3,.sinus
 shufps   xmm3,xmm3,01010101b ; xm3 = brdcs cos
 subps	  xmm4,xmm3	      ; xm4 = brdcst 1 - cos
 shufps   xmm5,xmm5,01000000b ; xm5 = l->h: nxxxy
 shufps   xmm6,xmm6,01011001b ; xm6 = l->h: nyzyy

 movhlps  xmm7,xmm4
 movlhps  xmm7,xmm0	      ; xm7 = lo -> hi: 1-cos, 1-cos, 1.0, -1.0
 mulps	  xmm6,xmm5
 shufps   xmm7,xmm7,10000000b ; xm7 = lo -> hi: 1-c 1-c 1-c 1.0
 mulps	  xmm6,xmm7	      ; xm6 = 1st part of add/sub 2, 3 ,4, 5 mx position
 movlps   xmm5,.ny
 movlhps  xmm5,xmm2	      ; xm5 = lo -> hi: ny, nz, 1-nx^2, 1-ny^2
 movlps   xmm7,.sinus
 shufps   xmm5,xmm5,11010001b ; xm5 = lo -> hi: nz, ny, nz, 1-ny^2

 shufps   xmm7,xmm7,01000000b ; xm7 = lo -> hi: s s s c
 mulps	  xmm7,xmm5	      ; xm7 = sec part of add/sub 2345 pos
 shufps   xmm0,xmm0,00010100b ; xm0= l->h: 1.0, -1.0, -1.0, 1.0
 mulps	  xmm7,xmm0
 addps	  xmm7,xmm6


 mulps	  xmm3,xmm2
 addps	  xmm3,xmm1

 movlps   [edi],xmm3
			      ; all from 2345
			      ; seems be ok
 movups   [edi+4],xmm7
 movups   xmm5,.nx
 movhlps  xmm7,xmm4
 movlhps  xmm7,xmm0
 movaps   xmm6,xmm5
 shufps   xmm5,xmm5,10101010b ; brdcs nz
 shufps   xmm6,xmm6,10010001b ; l->h: nyxyz
 shufps   xmm7,xmm7,10000000b ; l->h: 1-cos, 1-cos, 1-cos, 1.0
 mulps	  xmm5,xmm6
 mulps	  xmm5,xmm7	      ; 1st art of add/sub 6, 7, 8, 9 mx pos
 movhps   xmm6,.nx
 movhlps  xmm6,xmm2	      ; xm2 = 1 - nxyz?^2
 movlps   xmm7,.sinus
 shufps   xmm6,xmm6,00101110b ; xm6 = l->h: nx, ny, nx, 1-nz^2
 shufps   xmm7,xmm7,01000000b
 shufps   xmm0,xmm0,00100000b
 mulps	  xmm6,xmm7
 mulps	  xmm6,xmm0
 addps	  xmm6,xmm5
 movups   [edi+20],xmm6

 add	  esp,8

ret
